package cn.virens.dd.utils;

import java.util.Date;

import cn.hutool.core.date.DateUtil;

public class SunTimesUtil {
	private final static double h = -0.833;//日出日落时太阳的位置
	private final static double UTo = 180.0;//上次计算的日落日出时间，初始迭代值180.0

	private final static long hour_ms = 60 * 60 * 1000;
	private final static long days_ms = 24 * 60 * 60 * 1000;
	private final static Date d2000 = DateUtil.parseDate("2000-01-01");

	//求从格林威治时间公元2000年1月1日到计算日天数dayss
	private static long days(Date date) {
		return (date.getTime() - d2000.getTime()) / days_ms;
	}

	//求格林威治时间公元2000年1月1日到计算日的世纪数t
	private static double t_century(long days, double UTo) {
		return ((double) days + UTo / 360) / 36525;
	}

	//求太阳的平黄径

	private static double L_sun(double t_century) {
		return (280.460 + 36000.770 * t_century);
	}

	//求太阳的平近点角
	private static double G_sun(double t_century) {
		return (357.528 + 35999.050 * t_century);
	}

	//求黄道经度
	private static double ecliptic_longitude(double L_sun, double G_sun) {
		return (L_sun + 1.915 * Math.sin(G_sun * Math.PI / 180) + 0.02 * Math.sin(2 * G_sun * Math.PI / 180));

	}

	//求地球倾角
	private static double earth_tilt(double t_century) {
		return (23.4393 - 0.0130 * t_century);
	}

	//求太阳偏差
	private static double sun_deviation(double earth_tilt, double ecliptic_longitude) {
		return (180 / Math.PI * Math.asin(Math.sin(Math.PI / 180 * earth_tilt) * Math.sin(Math.PI / 180 * ecliptic_longitude)));
	}

	//求格林威治时间的太阳时间角GHA
	private static double GHA(double UTo, double G_sun, double ecliptic_longitude) {
		return (UTo - 180 - 1.915 * Math.sin(G_sun * Math.PI / 180) - 0.02 * Math.sin(2 * G_sun * Math.PI / 180) + 2.466 * Math.sin(2 * ecliptic_longitude * Math.PI / 180) - 0.053 * Math.sin(4 * ecliptic_longitude * Math.PI / 180));

	}

	//求修正值e
	private static double e(double h, double glat, double sun_deviation) {
		return 180 / Math.PI * Math.acos((Math.sin(h * Math.PI / 180) - Math.sin(glat * Math.PI / 180) * Math.sin(sun_deviation * Math.PI / 180)) / (Math.cos(glat * Math.PI / 180) * Math.cos(sun_deviation * Math.PI / 180)));

	}

	//求日出时间
	private static double UT_rise(double UTo, double GHA, double glong, double e) {
		return (UTo - (GHA + glong + e));
	}

	//求日落时间
	private static double UT_set(double UTo, double GHA, double glong, double e) {
		return (UTo - (GHA + glong - e));
	}

	//判断并返回结果（日出）
	private static double result_rise(double UT, double UTo, double glong, double glat, long dayss) {
		if (Math.abs(UT - UTo) >= 0.1) {
			double t_century = t_century(dayss, UT);

			double G_sun = G_sun(t_century);
			double L_sun = L_sun(t_century);

			double earth_tilt = earth_tilt(t_century);
			double ecliptic_longitude = ecliptic_longitude(L_sun, G_sun);
			double sun_deviation = sun_deviation(earth_tilt, ecliptic_longitude);

			double e = e(h, glat, sun_deviation);
			double GHA = GHA(UT, G_sun, ecliptic_longitude);

			double _UT = UT_rise(UT, GHA, glong, e);

			return result_rise(_UT, UT, glong, glat, dayss);
		} else {
			return UT;
		}
	}

	//判断并返回结果（日落）
	private static double result_set(double UT, double UTo, double glong, double glat, long days) {
		if (Math.abs(UT - UTo) >= 0.1) {
			double t_century = t_century(days, UT);

			double G_sun = G_sun(t_century);
			double L_sun = L_sun(t_century);

			double earth_tilt = earth_tilt(t_century);
			double ecliptic_longitude = ecliptic_longitude(L_sun, G_sun);
			double sun_deviation = sun_deviation(earth_tilt, ecliptic_longitude);

			double e = e(h, glat, sun_deviation);
			double GHA = GHA(UT, G_sun, ecliptic_longitude);

			double _UT = UT_set(UT, GHA, glong, e);

			return result_set(_UT, UT, glong, glat, days);
		} else {
			return UT;
		}
	}

	/** 日出时刻 */
	public static Date getSunrise(double glong, double glat, Date sunTime) {
		long days = days(sunTime);
		double t_century = t_century(days, UTo);

		double G_sun = G_sun(t_century);
		double L_sun = L_sun(t_century);

		double earth_tilt = earth_tilt(t_century);
		double ecliptic_longitude = ecliptic_longitude(L_sun, G_sun);
		double sun_deviation = sun_deviation(earth_tilt, ecliptic_longitude);

		double e = e(h, glat, sun_deviation);
		double GHA = GHA(UTo, G_sun, ecliptic_longitude);

		double UT_rise = UT_rise(UTo, GHA, glong, e);
		double sunrise = result_rise(UT_rise, UTo, glong, glat, days);

		return toTime(sunrise);
	}

	/** 日落时刻 */
	public static Date getSunset(double glong, double glat, Date sunTime) {
		long days = days(sunTime);
		double t_century = t_century(days, UTo);

		double G_sun = G_sun(t_century);
		double L_sun = L_sun(t_century);

		double earth_tilt = earth_tilt(t_century);
		double ecliptic_longitude = ecliptic_longitude(L_sun, G_sun);
		double sun_deviation = sun_deviation(earth_tilt, ecliptic_longitude);

		double e = e(h, glat, sun_deviation);
		double GHA = GHA(UTo, G_sun, ecliptic_longitude);

		double UT_set = UT_set(UTo, GHA, glong, e);
		double sunset = result_set(UT_set, UTo, glong, glat, days);

		return toTime(sunset);
	}

	private static Date toTime(double value) {
		return new Date((long) (value / 15 * hour_ms));
	}

	public static void main(String[] args) {
		Date str1 = getSunrise(106.5234456135, 29.325124, DateUtil.parse("2020-09-30"));
		Date str2 = getSunset(106.5234456135, 29.325124, DateUtil.parse("2020-09-30"));

		System.out.println("日出时间：" + DateUtil.formatDateTime(str1));
		System.out.println("日落时间：" + DateUtil.formatDateTime(str2));
	}
}